#include "spline.h"

//X e Y son los puntos a interpolar respectivamente
//en matCoef se devuelve el resultado
void spline(const vector<double>& X,const vector<double>& Y,vector< vector<double> >& matCoef)
{
	// Spline natural
	// S(X)j = aj + b(x - x0) + c(x  - x0)^2 +  d(x - x0)^3

	int n;
	int i, j;
	vector<double> h, alpha, u, z, l, c;
	
	//assert(X.size() == Y.size() && (int) X.size() >= 2);

	n = (int) X.size() - 1;
	h.resize(n);

	//PASO 1
	for (i = 0; i <= n - 1; i++)
		h[i] = X[i + 1] - X[i];

	alpha.resize(n + 1);
	u.resize(n);
	z.resize(n + 1);
	l.resize(n + 1);
	c.resize(n + 1);

	alpha[0] = 0;
	alpha[n] = 0;

	//PASO 2
	for (i = 1; i <= n - 1; i++)
		alpha[i] = (3.0 / h[i]) * (Y[i + 1] - Y[i]) - (3.0 / h[i - 1]) * (Y[i] - Y[i - 1]);

	//PASO 3
	l[0] = 1;
	u[0] = 0;
	z[0] = 0;

	//PASO 4
	for (i = 1; i <= n - 1; i++)
	{
		l[i] = 2.0 * (X[i + 1] - X[i - 1]) - h[i - 1] * u[i - 1];
		u[i] = h[i] / l[i];
		z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i];
	}

	//PASO 5
	l[n] = 1;
	z[n] = 0;
	c[n] = 0; 

	matCoef.resize(n);
	for (i = 0; i < (int) matCoef.size(); i++)
		matCoef[i].resize(4);

	//PASO 6
	for (j = n - 1; j >= 0; j--)
	{
		matCoef[j][0] = Y[j];  // A
		
		c[j] = z[j] - u[j] * c[j + 1]; // C
		matCoef[j][1] = (Y[j + 1] - Y[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3.0; // B
		matCoef[j][3] = (c[j + 1] - c[j]) / (3.0 * h[j]); // D
	}

	for (j = n - 1; j >= 0; j--)
		matCoef[j][2] = c[j]; // C
}

double calcularValor(vector< vector<double> >&miSpline,int diff, int pos)
{
	double elA=miSpline[pos][0];
	double elB=miSpline[pos][1]*(diff);
	double elC=miSpline[pos][2]*(pow(diff,2));
	double elD=miSpline[pos][3]*(pow(diff,3));
	return elA+elB+elC+elD;
}
